Давайте уточним правило трёх сигм. Утверждение: 99.7% вероятностной массы случайной величины X∼N(μ,σ2) лежит в интервале μ±c⋅σ. Чему равно точное значение константы c? Округлите ответ до четырёх знаков после десятичной точки.
В пятилетнем рандомизированном исследовании Гарвардской медицинской школы 11037 испытуемых через день принимали аспирин, а ещё 11034 — плацебо. Исследование было слепым, то есть, испытуемые не знали, что именно они принимают. За 5 лет инфаркт случился у 104 испытуемых, принимавших аспирин, и у 189 принимавших плацебо. Оцените, насколько вероятность инфаркта снижается при приёме аспирина. Округлите ответ до четырёх знаков после десятичной точки.
Постройте теперь 95% доверительный интервал для снижения вероятности инфаркта при приёме аспирина. Чему равна его верхняя граница? Округлите ответ до четырёх знаков после десятичной точки.
Продолжим анализировать данные эксперимента Гарвардской медицинской школы. Для бернуллиевских случайных величин X∼Ber(p) часто вычисляют величину p1−p, которая называется шансами (odds). Чтобы оценить шансы по выборке, вместо p нужно подставить p^. Например, шансы инфаркта в контрольной группе, принимавшей плацебо, можно оценить как ≈0.0174 Оцените, во сколько раз понижаются шансы инфаркта при регулярном приёме аспирина. Округлите ответ до четырёх знаков после десятичной точки.
Величина, которую вы оценили в предыдущем вопросе, называется отношением шансов. Постройте для отношения шансов 95% доверительный интервал с помощью бутстрепа. Чему равна его нижняя граница? Округлите ответ до 4 знаков после десятичной точки. Чтобы получить в точности такой же доверительный интервал, как у нас: составьте векторы исходов в контрольной и тестовой выборках так, чтобы в начале шли все единицы, а потом все нули; установите random seed=0; сделайте по 1000 псевдовыборок из каждой группы пациентов.
In [1]:
import scipy.stats
scipy.stats.norm.ppf(0.9985)
Out[1]:
In [12]:
import numpy as np
n1=11037
k1=104
p1=float(k1)/n1
n2=11034
k2=189
p2=float(k2)/n2
p2-p1+1.96*np.sqrt(p2*(1-p2)/n2+p1*(1-p1)/n1)
Out[12]:
In [6]:
odds1=p1/(1-p1)
odds2=p2/(1-p2)
print odds1, odds2
print odds2/odds1
In [16]:
def get_bootstrap_samples(data, n_samples):
indices = np.random.randint(0, len(data), (n_samples, len(data)))
samples = data[indices]
return samples
In [17]:
def stat_intervals(stat, alpha):
boundaries = np.percentile(stat, [100 * alpha / 2., 100 * (1 - alpha / 2.)])
return boundaries
In [18]:
def get_p(k,n):
return float(k)/n
In [19]:
def get_odds(p):
return p/(1-p)
In [31]:
def get_chances_relation(sample1, sample2):
return get_odds(get_p(sum(sample2), len(sample2)))/get_odds(get_p(sum(sample1), len(sample1)))
In [32]:
data1 = np.zeros(n1)
data1[:k1]=1
data2 = np.zeros(n2)
data2[:k2]=1
In [33]:
np.random.seed(0)
samples1 = get_bootstrap_samples(data1, 1000)
samples2 = get_bootstrap_samples(data2, 1000)
rel_scores = []
for i in xrange(1000):
rel_scores.append(get_chances_relation(samples1[i], samples2[i]))
In [34]:
print "95% confidence interval:", stat_intervals(rel_scores, 0.05)
In [ ]: